**************************************************************
* Simulation with actual MSCZ and BCP data
**************************************************************
clear all
* RUN 00_path_master.do FIRST TO GET FILEPATHS

set seed 379215

**************************************************************
* 1. Define global parameters
**************************************************************

local simulations = 500
local leads = 3
local lags = 3
local total_coefs = `leads' + `lags' + 1  // 3 leads, 3 lags, 1 contemporary
local long_total_coefs = `total_coefs' + 8 // 1990 - 1997

frame create msczes
frame change msczes

use "$intermediate/qcew_msczes.dta", clear

frame create bcpes
frame change bcpes

use "$intermediate/qcew_bcpes.dta", clear

* Create matrices to store results
matrix twfe_static_mscz = J(`simulations', 1, .)
matrix twfe_static_mscz_00 = J(`simulations', 1, .)
matrix twfe_static_bcp = J(`simulations', 1, .)
matrix twfe_static_bcp_00 = J(`simulations', 1, .)
matrix twfe_dl_mscz = J(`simulations', `total_coefs', .)
matrix twfe_dl_bcp = J(`simulations', `total_coefs', .)
matrix xtevent_mscz = J(`simulations', `total_coefs', .)
matrix xtevent_bcp = J(`simulations', `total_coefs', .)
matrix lpdidlong_mscz = J(`simulations', `long_total_coefs', .)
matrix lpdidlong_bcp = J(`simulations', `long_total_coefs', .)
matrix lpdidshort_mscz = J(`simulations', `total_coefs', .)
matrix lpdidshort_bcp = J(`simulations', `total_coefs', .)


**************************************************************
* MSCZ
**************************************************************

frame change msczes


g post2000 = year>=2000 
*g post2000 = year > 2000
g post1997 = year>1997
gegen event_post2000 = max(event_combined*post2000), by(pair_id_czstate)


g trend = 0
replace trend = 0.087  if year == 1990 & event_post2000 == 1
replace trend = 0.086  if year == 1991 & event_post2000 == 1
replace trend = 0.101  if year == 1992 & event_post2000 == 1
replace trend = 0.08  if year == 1993 & event_post2000 == 1
replace trend = 0.085  if year == 1994 & event_post2000 == 1
replace trend = 0.081  if year == 1995 & event_post2000 == 1
replace trend = 0.089  if year == 1996 & event_post2000 == 1
replace trend = 0.058  if year == 1997 & event_post2000 == 1

xtset pair_id_czstate year

sort pair_id_czstate year, stable
gen lead1 = D.F1.lmw
gen lead2 = D.F2.lmw
gen lead3 = D.F3.lmw
forval j=1/3 {
	replace lead`j' = 0 if lead`j' == .
}

gen lag1  = D.L1.lmw
gen lag2  = D.L2.lmw
gen lag3  = L3.lmw
forval j=1/3 {
	replace lag`j' = 0 if lag`j' == .
}
	
g dlmw = D.lmw
replace dlmw=0 if dlmw==.

gen mw_rise = mw > L.mw

gegen pair_id_year = group(pair_id_num year)


* estimate SD for error *
reghdfe   lepop75  if inrange(year,1990,2016), absorb(pair_id_czstate pair_id_num##year) cluster(statefips) res
predict yres, res
sum yres
scalar errorSD=r(sd)
drop yres


forval sim = 1/`simulations' {
	
	/****** Draw Outcome ******/
	g unitFE = rnormal(0,0.03)
	bys pair_id_czstate (year): replace unitFE = unitFE[1]

	g timeFE = rnormal(0, 0.03)
	bys pair_id_year (pair_id_czstate): replace timeFE = timeFE[1]
	
	g error = rnormal(0,errorSD)
 
	g trend_error = 0
	replace trend_error= rnormal(0,0.015) if year<=1997 & event_post2000==1
	bys event_post2000 post1997 (pair_id_czstate year): replace trend_error = trend_error[1]

	g outcome = trend +  trend_error +  timeFE   + unitFE  +  error 
	
	/***************************/
	
	
	sort pair_id_czstate year, stable
	forvalues i = 1/3 {
		gen outcome_pre`i' = L`i'.outcome
		gen outcome_post`i' = F`i'.outcome
	}
	
	loc j = 1
	
	forvalues i = 3(-1)1 {
		gen outcome_dif`j' = outcome_pre`i' - outcome_pre1
		loc ++j
	}
	
	gen outcome_dif4 = outcome - outcome_pre1
	
	loc j = 5
	
	forvalues i = 1/3{
		gen outcome_dif`j' = outcome_post`i' - outcome_pre1
		loc ++j
	}
	
	forvalues i = 1/7 {
		local j = `i' + 8
		rename outcome_dif`i' outcome_dif`j'
	}
	
	forvalues i = 1990/1997 {
			
			gen _outcome_longpre = outcome if year == `i'
			bys pair_id_czstate (year): egen outcome_longpre`i' = mean(_outcome_longpre)
			drop _outcome_longpre
	}
	
	forvalues i = 1/8 {
		local k = `i' + 1989
		gen outcome_dif`i' = outcome_longpre`k' - outcome_pre1
	}
	
	forvalues i = 1/15 {
				replace outcome_dif`i' = . if inlist(year,1990,1991,1996,1997,2007,2008,2009)
	}
			
	replace outcome_dif15 = . if inrange(year,1993,1995) | inrange(year,2004,2006) | inrange(year,2017,2019)
	replace outcome_dif14 = . if inrange(year,1994,1995) | inrange(year,2005,2006) | inrange(year,2018,2019)
	replace outcome_dif13 = . if inlist(year,1995,2006,2019)
	replace outcome_dif8 = . if year == 2000	
		
		
		
	reghdfe outcome lead2 lead1 dlmw lag1 lag2 lag3 if inrange(year,1990,2016), absorb(pair_id_czstate pair_id_num##year) cluster(statefips)
	matrix twfe_dl_mscz[`sim', 1] = -_b[lead1]
	matrix twfe_dl_mscz[`sim', 2] = _b[lead2] - _b[lead1]
	matrix twfe_dl_mscz[`sim', 3] = _b[lead1] - _b[lead1]
	matrix twfe_dl_mscz[`sim', 4] = _b[dlmw] - _b[lead1]
	matrix twfe_dl_mscz[`sim', 5] = _b[lag1] - _b[lead1]
	matrix twfe_dl_mscz[`sim', 6] = _b[lag2] - _b[lead1]
	matrix twfe_dl_mscz[`sim', 7] = _b[lag3] - _b[lead1]
		
	reghdfe outcome lmw if inrange(year,1990,2016), absorb(pair_id_czstate pair_id_num##year) cluster(statefips)
	matrix twfe_static_mscz[`sim', 1] = _b[lmw]
	
	reghdfe outcome lmw if inrange(year,2000,2016), absorb(pair_id_czstate pair_id_num##year) cluster(statefips)
	matrix twfe_static_mscz_00[`sim', 1] = _b[lmw]

	forvalues i = 1/15 {
		reghdfe outcome_dif`i' event_combined if (event_combined == 1 | clean_control == 1) & inrange(year,2000,2016), a(pair_id_num##year) cluster(statefips)
		matrix lpdidlong_mscz[`sim', `i'] = _b[event_combined]
	}
	
	forvalues i = 1/7 {
		local j = `i' + 8
		reghdfe outcome_dif`j' event_combined if (event_combined == 1 | clean_control == 1) & inrange(year,1990,2016), a(year) cluster(statefips)
		matrix lpdidshort_mscz[`sim', `i'] = _b[event_combined]
	}
		
		xtevent outcome if inrange(year,1990,2016), w(2)  pol(mw_rise) nofe  ///
					note cluster(statefips) reghdfe  addabsorb(pair_id_czstate pair_id_num##year) /*impute(stag)*/
		matrix xtevent_mscz[`sim', 1] = _b[_k_eq_m3]
        matrix xtevent_mscz[`sim', 2] = _b[_k_eq_m2]
        matrix xtevent_mscz[`sim', 3] = 0
        matrix xtevent_mscz[`sim', 4] = _b[_k_eq_p0]
        matrix xtevent_mscz[`sim', 5] = _b[_k_eq_p1]
        matrix xtevent_mscz[`sim', 6] = _b[_k_eq_p2]
        matrix xtevent_mscz[`sim', 7] = _b[_k_eq_p3]
		
		drop unitFE timeFE error outcome* trend_error
    }

    **************************************************************
    * Convert results to data sets and define event time
    **************************************************************
    svmat twfe_dl_mscz, names(twfe_dl)
    svmat twfe_static_mscz, names(twfe_static)
	svmat twfe_static_mscz_00, names(twfe_static_00)
	svmat lpdidlong_mscz, names(lplong)
	svmat lpdidshort_mscz, names(lpshort)
    svmat xtevent_mscz, names(xtevent)

    gen lag = _n - 4

    **************************************************************
    * Compute means and SDs for each estimator
    **************************************************************
    * TWFE
    gen coef_mean_twfe = .
    gen coef_sd_twfe   = .
    forval i = 1/7 {
        summarize twfe_dl`i', detail
        replace coef_mean_twfe = r(mean) in `i'
        replace coef_sd_twfe   = r(sd)   in `i'
    }
    gen twfe_lower = coef_mean_twfe - 1.96 * coef_sd_twfe
    gen twfe_upper = coef_mean_twfe + 1.96 * coef_sd_twfe
	
	 * TWFE - static
    gen coef_mean_twfestatic = .
    gen coef_sd_twfestatic   = .
    
	summarize twfe_static1, detail
	replace coef_mean_twfestatic = r(mean) in 1
	replace coef_sd_twfestatic   = r(sd)   in 1
	scalar twfestatic_mean = int(100*r(mean))/100
	scalar twfestatic_SD = int(100*r(sd))/100
	
	gen coef_mean_twfestatic_00 = .
    gen coef_sd_twfestatic_00   = .
    
	summarize twfe_static_001, detail
	replace coef_mean_twfestatic_00 = r(mean) in 1
	replace coef_sd_twfestatic_00   = r(sd)   in 1
	scalar twfestatic_00_mean = int(100*r(mean))/100
	scalar twfestatic_00_SD = int(100*r(sd))/100
	
	* LP-DiD long
    gen coef_mean_lpdidlong = .
    gen coef_sd_lpdidlong   = .
    forval i = 1/15 {
        summarize lplong`i', detail
        replace coef_mean_lpdidlong = r(mean) in `i'
        replace coef_sd_lpdidlong   = r(sd)   in `i'
    }
    gen lplong_lower = coef_mean_lpdidlong - 1.96 * coef_sd_lpdidlong
    gen lplong_upper = coef_mean_lpdidlong + 1.96 * coef_sd_lpdidlong

	gen coef_mean_lpdidshort = .
    gen coef_sd_lpdidshort   = .
    forval i = 1/7 {
        summarize lpshort`i', detail
        replace coef_mean_lpdidshort = r(mean) in `i'
        replace coef_sd_lpdidshort = r(sd)   in `i'
    }
    gen lpshort_lower = coef_mean_lpdidshort - 1.96 * coef_sd_lpdidshort
    gen lpshort_upper = coef_mean_lpdidshort + 1.96 * coef_sd_lpdidshort
	
	/*
	* LP-DiD short *
	
	cap drop *short*
	
	g coef_mean_lpdidshort = coef_mean_lpdidlong[_n+8]
	g coef_sd_lpdidshort = coef_sd_lpdidlong[_n+8]
	g lpshort_lower = lplong_lower[_n+8]
	g lpshort_upper = lplong_upper[_n+8]
*/
	* XTEVENT
    gen coef_mean_xtevent = .
    gen coef_sd_xtevent   = .
    forval i = 1/7 {
        summarize xtevent`i', detail
        replace coef_mean_xtevent = r(mean) in `i'
        replace coef_sd_xtevent   = r(sd)   in `i'
    }
    gen xtevent_lower = coef_mean_xtevent - 1.96 * coef_sd_xtevent
    gen xtevent_upper = coef_mean_xtevent + 1.96 * coef_sd_xtevent
	
	
	**************************************************************
    * TWFE-logMW graph
    **************************************************************
	
        twoway (line coef_mean_twfe lag if lag <=3 & lag>=-3, lcolor(blue%50) lwidth(med)) ///
               (rarea twfe_lower twfe_upper lag if lag <=3 & lag>=-3, fcolor(blue%20) lcolor(blue%0)) ///
               (scatter coef_mean_twfe lag if lag <=3 & lag>=-3, mcolor(blue) msymbol(circle)), ///
			    xtitle("Event time", size(small)) ytitle("Cumulative response", size(small)) ///
               legend(off) xlab(-3(1)3) yline(0,lcolor(gs8) lpattern(dash)) ///
               ysize(3) xsize(4) name(sim_twfedl_mscz, replace) ylab(-.25(.05).2) ///
			   text( -.2 0 "Static TWFE est (SD): `=twfestatic_mean'(`=twfestatic_SD')" ) title("TWFE with log MW")

	**************************************************************
    * LP-DiD and XTEVENT graph
    **************************************************************

			   
		twoway (line coef_mean_lpdidshort lag if lag <=3 & lag>=-3, lcolor(green%50) lwidth(med)) ///
			   (scatter coef_mean_lpdidshort lag if lag <=3 & lag>=-3, mcolor(green%50) ) ///
               (rarea lpshort_lower lpshort_upper lag if lag <=3 & lag>=-3, fcolor(green%20) lcolor(green%0)) ///
			   (rarea xtevent_lower xtevent_upper lag if lag <=3 & lag>=-3, fcolor(orange%20) lcolor(orange%20)) ///
               (scatter coef_mean_xtevent lag if lag <=3 & lag>=-3, mcolor(orange) msymbol(Dh))  ///
			   (line coef_mean_xtevent lag if lag <=3 & lag>=-3, lcolor(orange) ), ///
			   xtitle("Event time", size(small)) ytitle("Cumulative response", size(small)) ///
			   yline(0,lcolor(gs8) lpattern(dash)) ///
               legend(order(2 "LP-DiD (Combined state-events)" 5 "XTEVENT (All changes)") row(2) pos(6)) xlab(-3(1)3) ///
               ysize(3) xsize(4) name(sim_eventstudy_mscz, replace) ylab(-.04(.01).02) title("LP-DiD event study (all combined state events)," "and XTEVENT (binary changes)")
			   
	
 	
	**************************************************************
    * Extended pre-trends graph
    **************************************************************
	
	insobs 1, a(8)
	local last_90s = coef_mean_lpdidlong[8]
	local first_event = coef_mean_lpdidlong[10]
	replace lag = _n - 4
	
	twoway (line coef_mean_lpdidlong lag if lag>=6 & lag<=12, lcolor(green*0.5) lpattern(solid)) ///
	(rcap lplong_lower lplong_upper lag if lag>=6 & lag<=12, lcolor(green*0.5) lpattern(solid)) ///
	(line coef_mean_lpdidlong lag if lag<5, lcolor(green*0.25) lpattern(solid)) ///
	(rcap lplong_lower lplong_upper lag if lag<5, lcolor(green*0.25) lpattern(solid)) ///
	(scatteri `last_90s' 4 `first_event' 6, recast(line) lcolor(green*0.25) lpattern(shortdash) lwidth(med)) ///
	, xlabel(-3 "1990" -1 "1992" 1 "1994" 3 "1996" 6 "-3" 8 "-1" 10 "1" 12 "3") ///
	text(-0.085 1 "Years in 1990s", color(green*0.25) size(small)) ///
	text(-0.095 1 "(Extended pre-period)", color(green*0.25) size(small)) ///
	text(-0.085 9 "Years around MW increase event", color(green*0.5) size(small)) ///
	text(-0.095 9 "(Standard event window)", color(green*0.5) size(small)) ///
	xtitle(" ") ylab(-0.04(0.04)0.12) yline(0,lcolor(gs8) lpattern(dash)) ///
	ytitle("Change in outcome", size(small)) legend(off) name(sim_pret_mscz, replace)  title("LP-DiD event study with extended" "pre-period (post-2000 combined state events)")

	
	
 	 
	
**************************************************************
* BCP
**************************************************************

frame change bcpes


g post2000 = year>=2000 
g post1997 = year>1997
gegen event_post2000 = max(event_combined*post2000), by(pair_id_county)


g trend = 0
replace trend =  0.031   if year == 1990 & event_post2000 == 1
replace trend =  0.033  if year == 1991 & event_post2000 == 1
replace trend =  0.024  if year == 1992 & event_post2000 == 1
replace trend =  0.017  if year == 1993 & event_post2000 == 1
replace trend =  0.021  if year == 1994 & event_post2000 == 1
replace trend =  0.013  if year == 1995 & event_post2000 == 1
replace trend =  0.007  if year == 1996 & event_post2000 == 1
replace trend =  -0.003  if year == 1997 & event_post2000 == 1


xtset pair_id_county year

sort pair_id_county year, stable
gen lead1 = D.F1.lmw
gen lead2 = D.F2.lmw
gen lead3 = D.F3.lmw
forval j=1/3 {
	replace lead`j' = 0 if lead`j' == .
}

gen lag1  = D.L1.lmw
gen lag2  = D.L2.lmw
gen lag3  = L3.lmw
forval j=1/3 {
	replace lag`j' = 0 if lag`j' == .
}
	
g dlmw = D.lmw
replace dlmw=0 if dlmw==.

gen mw_rise = mw > L.mw

gegen pair_id_year = group(pair_id_num year)



* estimate SD for error *
reghdfe   lepop75  if inrange(year,1990,2016), absorb(pair_id_county pair_id_num##year) cluster(statefips) res
predict yres, res
sum yres
scalar errorSD=r(sd)
drop yres


forval sim = 1/`simulations' {
	g unitFE = rnormal(0,0.03)
	bys pair_id_county (year): replace unitFE = unitFE[1]

	g timeFE = rnormal(0, 0.03)
	bys pair_id_year (pair_id_county): replace timeFE = timeFE[1]
	
	g error = rnormal(0,errorSD)
 
	
	g trend_error = 0
	replace trend_error= rnormal(0,0.015) if year<=1997 & event_post2000==1
	bys event_post2000 post1997 (pair_id_county year): replace trend_error = trend_error[1]

	g outcome = trend +  trend_error +  timeFE   + unitFE  +  error 
	
	 
	sort pair_id_county year, stable
	forvalues i = 1/3 {
		gen outcome_pre`i' = L`i'.outcome
		gen outcome_post`i' = F`i'.outcome
	}
	
	loc j = 1
	
	forvalues i = 3(-1)1 {
		gen outcome_dif`j' = outcome_pre`i' - outcome_pre1
		loc ++j
	}
	
	gen outcome_dif4 = outcome - outcome_pre1
	
	loc j = 5
	
	forvalues i = 1/3{
		gen outcome_dif`j' = outcome_post`i' - outcome_pre1
		loc ++j
	}
	
	forvalues i = 1/7 {
		local j = `i' + 8
		rename outcome_dif`i' outcome_dif`j'
	}
	
	forvalues i = 1990/1997 {
			
			gen _outcome_longpre = outcome if year == `i'
			bys pair_id_county (year): egen outcome_longpre`i' = mean(_outcome_longpre)
			drop _outcome_longpre
	}
	
	forvalues i = 1/8 {
		local k = `i' + 1989
		gen outcome_dif`i' = outcome_longpre`k' - outcome_pre1
	}
	
	forvalues i = 1/15 {
				replace outcome_dif`i' = . if inlist(year,1990,1991,1996,1997,2007,2008,2009)
	}
			
	replace outcome_dif15 = . if inrange(year,1993,1995) | inrange(year,2004,2006) | inrange(year,2017,2019)
	replace outcome_dif14 = . if inrange(year,1994,1995) | inrange(year,2005,2006) | inrange(year,2018,2019)
	replace outcome_dif13 = . if inlist(year,1995,2006,2019)
	replace outcome_dif8 = . if year == 2000	
		
		
		
	reghdfe outcome lead2 lead1 dlmw lag1 lag2 lag3 if inrange(year,1990,2016), absorb(pair_id_county pair_id_num##year) cluster(statefips)
	matrix twfe_dl_bcp[`sim', 1] = -_b[lead1]
	matrix twfe_dl_bcp[`sim', 2] = _b[lead2] - _b[lead1]
	matrix twfe_dl_bcp[`sim', 3] = _b[lead1] - _b[lead1]
	matrix twfe_dl_bcp[`sim', 4] = _b[dlmw] - _b[lead1]
	matrix twfe_dl_bcp[`sim', 5] = _b[lag1] - _b[lead1]
	matrix twfe_dl_bcp[`sim', 6] = _b[lag2] - _b[lead1]
	matrix twfe_dl_bcp[`sim', 7] = _b[lag3] - _b[lead1]
		
	reghdfe outcome lmw if inrange(year,1990,2016), absorb(pair_id_county pair_id_num##year) cluster(statefips)
	matrix twfe_static_bcp[`sim', 1] = _b[lmw]
	
	reghdfe outcome lmw if inrange(year,2000,2016), absorb(pair_id_county pair_id_num##year) cluster(statefips)
	matrix twfe_static_bcp_00[`sim', 1] = _b[lmw]

	forvalues i = 1/15 {
		reghdfe outcome_dif`i' event_combined if (event_combined == 1 | clean_control == 1) & inrange(year,2000,2016), a(pair_id_num##year) cluster(statefips)
		matrix lpdidlong_bcp[`sim', `i'] = _b[event_combined]
	}
	
	forvalues i = 1/7 {
		local j = `i' + 8
		reghdfe outcome_dif`j' event_combined if (event_combined == 1 | clean_control == 1) & inrange(year,1990,2016), a(year) cluster(statefips)
		matrix lpdidshort_bcp[`sim', `i'] = _b[event_combined]
	}
		
		xtevent outcome if inrange(year,1990,2016) , w(2)  pol(mw_rise) nofe  ///
					note cluster(statefips) reghdfe  addabsorb(pair_id_county pair_id_num##year) /*impute(stag)*/
		matrix xtevent_bcp[`sim', 1] = _b[_k_eq_m3]
        matrix xtevent_bcp[`sim', 2] = _b[_k_eq_m2]
        matrix xtevent_bcp[`sim', 3] = 0
        matrix xtevent_bcp[`sim', 4] = _b[_k_eq_p0]
        matrix xtevent_bcp[`sim', 5] = _b[_k_eq_p1]
        matrix xtevent_bcp[`sim', 6] = _b[_k_eq_p2]
        matrix xtevent_bcp[`sim', 7] = _b[_k_eq_p3]
		
		drop unitFE timeFE error trend_error outcome*
    }

    **************************************************************
    * Convert results to data sets and define event time
    **************************************************************
    svmat twfe_dl_bcp, names(twfe_dl)
    svmat twfe_static_bcp, names(twfe_static)
	svmat twfe_static_bcp_00, names(twfe_static_00)
	svmat lpdidlong_bcp, names(lplong)
	svmat lpdidshort_bcp, names(lpshort)
    svmat xtevent_bcp, names(xtevent)

    gen lag = _n - 4

    **************************************************************
    * Compute means and SDs for each estimator
    **************************************************************
    * TWFE
    gen coef_mean_twfe = .
    gen coef_sd_twfe   = .
    forval i = 1/7 {
        summarize twfe_dl`i', detail
        replace coef_mean_twfe = r(mean) in `i'
        replace coef_sd_twfe   = r(sd)   in `i'
    }
    gen twfe_lower = coef_mean_twfe - 1.96 * coef_sd_twfe
    gen twfe_upper = coef_mean_twfe + 1.96 * coef_sd_twfe
	
	 * TWFE - static
    gen coef_mean_twfestatic = .
    gen coef_sd_twfestatic   = .
    
	summarize twfe_static1, detail
	replace coef_mean_twfestatic = r(mean) in 1
	replace coef_sd_twfestatic   = r(sd)   in 1
	scalar twfestatic_mean = int(100*r(mean))/100
	scalar twfestatic_SD = int(100*r(sd))/100
	
	gen coef_mean_twfestatic_00 = .
    gen coef_sd_twfestatic_00   = .
    
	summarize twfe_static_001, detail
	replace coef_mean_twfestatic_00 = r(mean) in 1
	replace coef_sd_twfestatic_00   = r(sd)   in 1
	scalar twfestatic_00_mean = int(100*r(mean))/100
	scalar twfestatic_00_SD = int(100*r(sd))/100
	
	    * LP-DiD long
    gen coef_mean_lpdidlong = .
    gen coef_sd_lpdidlong   = .
    forval i = 1/15 {
        summarize lplong`i', detail
        replace coef_mean_lpdidlong = r(mean) in `i'
        replace coef_sd_lpdidlong   = r(sd)   in `i'
    }
    gen lplong_lower = coef_mean_lpdidlong - 1.96 * coef_sd_lpdidlong
    gen lplong_upper = coef_mean_lpdidlong + 1.96 * coef_sd_lpdidlong
	
	gen coef_mean_lpdidshort = .
    gen coef_sd_lpdidshort   = .
    forval i = 1/7 {
        summarize lpshort`i', detail
        replace coef_mean_lpdidshort = r(mean) in `i'
        replace coef_sd_lpdidshort = r(sd)   in `i'
    }
    gen lpshort_lower = coef_mean_lpdidshort - 1.96 * coef_sd_lpdidshort
    gen lpshort_upper = coef_mean_lpdidshort + 1.96 * coef_sd_lpdidshort
	
/*	* LP-DiD short *
	
	cap drop *short*
	
	g coef_mean_lpdidshort = coef_mean_lpdidlong[_n+8]
	g coef_sd_lpdidshort = coef_sd_lpdidlong[_n+8]
	g lpshort_lower = lplong_lower[_n+8]
	g lpshort_upper = lplong_upper[_n+8]
*/

	* XTEVENT
    gen coef_mean_xtevent = .
    gen coef_sd_xtevent   = .
    forval i = 1/7 {
        summarize xtevent`i', detail
        replace coef_mean_xtevent = r(mean) in `i'
        replace coef_sd_xtevent   = r(sd)   in `i'
    }
    gen xtevent_lower = coef_mean_xtevent - 1.96 * coef_sd_xtevent
    gen xtevent_upper = coef_mean_xtevent + 1.96 * coef_sd_xtevent
	
	
	
	**************************************************************
    * TWFE-logMW graph
    **************************************************************
	
        twoway (line coef_mean_twfe lag if lag <=3 & lag>=-3, lcolor(blue%50) lwidth(med)) ///
               (rarea twfe_lower twfe_upper lag if lag <=3 & lag>=-3, fcolor(blue%20) lcolor(blue%0)) ///
               (scatter coef_mean_twfe lag if lag <=3 & lag>=-3, mcolor(blue) msymbol(circle)), ///
			    xtitle("Event time", size(small)) ytitle("Cumulative response", size(small)) ///
               legend(off) xlab(-3(1)3) yline(0,lcolor(gs8) lpattern(dash)) ///
               ysize(3) xsize(4) name(sim_twfedl_bcp, replace) ylab(-.25(.05).2) ///
			   text( -.2 0 "Static TWFE est (SD): `=twfestatic_mean'(`=twfestatic_SD')" ) title("TWFE with log MW")

	**************************************************************
    * LP-DiD and XTEVENT graph
    **************************************************************

			   
		twoway (line coef_mean_lpdidshort lag if lag <=3 & lag>=-3, lcolor(green%50) lwidth(med)) ///
			   (scatter coef_mean_lpdidshort lag if lag <=3 & lag>=-3, mcolor(green%50) ) ///
               (rarea lpshort_lower lpshort_upper lag if lag <=3 & lag>=-3, fcolor(green%20) lcolor(green%0)) ///
			   (rarea xtevent_lower xtevent_upper lag if lag <=3 & lag>=-3, fcolor(orange%20) lcolor(orange%20)) ///
               (scatter coef_mean_xtevent lag if lag <=3 & lag>=-3, mcolor(orange) msymbol(Dh))  ///
			   (line coef_mean_xtevent lag if lag <=3 & lag>=-3, lcolor(orange) ), ///
			    xtitle("Event time", size(small)) ytitle("Cumulative response", size(small)) ///
				yline(0,lcolor(gs8) lpattern(dash)) ///
               legend(order(2 "LP-DiD (Combined state-events)" 5 "XTEVENT (All changes)") row(2) pos(6)) xlab(-3(1)3) ///
               ysize(3) xsize(4) name(sim_eventstudy_bcp, replace) ylab(-.04(.01).02) title("LP-DiD event study (all combined state events)," "and XTEVENT (binary changes)")
			   
	
 	
	**************************************************************
    * Extended pre-trends graph
    **************************************************************
	
	insobs 1, a(8)
	local last_90s = coef_mean_lpdidlong[8]
	local first_event = coef_mean_lpdidlong[10]
	replace lag = _n - 4
	
	twoway (line coef_mean_lpdidlong lag if lag>=5 & lag<=12, lcolor(green*0.5) lpattern(solid)) ///
	(rcap lplong_lower lplong_upper lag if lag>=6 & lag<=12, lcolor(green*0.5) lpattern(solid)) ///
	(line coef_mean_lpdidlong lag if lag<5, lcolor(green*0.25) lpattern(solid)) ///
	(rcap lplong_lower lplong_upper lag if lag<5, lcolor(green*0.25) lpattern(solid)) ///
	(scatteri `last_90s' 4 `first_event' 6, recast(line) lcolor(green*0.25) lpattern(shortdash) lwidth(med)) ///
	, xlabel(-3 "1990" -1 "1992" 1 "1994" 3 "1996" 6 "-3" 8 "-1" 10 "1" 12 "3") ///
	text(-0.075 1 "Years in 1990s", color(green*0.25) size(small)) ///
	text(-0.085 1 "(Extended pre-period)", color(green*0.25) size(small)) ///
	text(-0.075 9 "Years around MW increase event", color(green*0.5) size(small)) ///
	text(-0.085 9 "(Standard event window)", color(green*0.5) size(small)) ///
	xtitle(" ")  ylab(-0.03(.03)0.15) yline(0,lcolor(gs8) lpattern(dash)) ///
	ytitle("Change in outcome", size(small)) legend(off) name(sim_pret_bcp, replace)  title("LP-DiD event study with extended" " pre-period (post-2000 combined state events)")
	

******* COMBINED GRAPH *********	
	 
 
graph combine 	sim_twfedl_mscz  sim_eventstudy_mscz sim_pret_mscz, col(1) name(MSCZ, replace) title(A. MSCZ)

graph combine 	sim_twfedl_bcp  sim_eventstudy_bcp sim_pret_bcp, col(1) name(BCP, replace) title(B. BCP)

gr combine MSCZ BCP, col(2) ysize(16) xsize(10)

graph export "$figures/Figure E1.pdf", replace	




	
 